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Abstract 

We consider the case when a set of spatially distributed sensors Qi,..., Q p make local observations, 
yy,... , y p , which are noisy versions of a signal of interest, x. Each sensor Qj transmits compressed 
information u ; about its measurements to the fusion center which should recover the original signal 
within a prescribed accuracy. Such an information processing relates to a wireless sensor network (WSN) 
scenario. The key problem is to find models of the sensors and fusion center so that they will be optimal 
in the sense of minimization of the associated error under a certain criterion, such as the mean square 
error (MSE). We determine the models from the technique which is a combination of the maximum 
block improvement (MBI) method JT], j2j and the generic Karhunen-Loeve transform (KLT) @ (based 
on the work in J4j, 0). Therefore, the proposed method unites the merits of both techniques fT), |[2l 
and a, ed, 0. As a result, our approach provides, in particular, the minimal MSE at each step of 
the version of the MBI method we use. The WSN model is represented in the form called the multi¬ 
compressor KLT-MBI transform. The multi-compressor KLT-MBI is given in terms of pseudo-inverse 
matrices and, therefore, it is numerically stable and always exists. In other words, the proposed WSN 
model provides compression, de-noising and reconstruction of distributed signals for the cases when 
known methods either are not applicable (because of singularity of associated matrices) or produce 
larger associated errors. Error analysis is provided. 
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Fig. 1. Block diagram of the WSN. Here, x is an estimation of x. 


Karhunen-Loeve transform, singular value decomposition, rank-reduced matrix approximation. 

I. Introduction 

A. Motivation 

We seek to find effective numerical algorithms for an information processing scenario that 
involves a set of spatially distributed sensors, Q 1} ..., Q p , and a fusion center, V. The sensors 
make local observations, yi,. ... y p , which are noisy versions of a signal of interest, x. Each 
sensor Q 3 transmits compressed information about its measurements to the fusion center 
which should recover the original signal within a prescribed accuracy. Such an information 
processing relates to a wireless sensor network (WSN) scenario. In the recent years, research 
and development on new and refined WSN techniques has increased at a remarkable rate (see, 
for example, J6]|, 0, 0, flU, iflOl . [HU . Ifl2l ). This is, in particular, because of a multitude of 
WSN applications due to their low deployment and maintenance cost. 

The key problem is to to find an effective way to compress and denoise each observation y 3 , 
where j = 1,,p, and then reconstruct all the compressed observations in the fusion center 
so that the reconstruction will be optimal in the sense of a minimization of the associated error 
under a certain criterion, such as the mean square error (MSE). A restriction is that the sensors 
cannot communicate with each other. Here, the term “compression” is treated in the same sense 
as in the known works on data compression (developed, for instance, in lfl3l . Ifl4t [[151 . |[T6l ). 
i.e. we say that observed signal y 3 with n 3 components is compressed if it is represented as 


August 20, 2015 


DRAFT 











3 


signal Uj with r 3 components where r 3 < nj, for j = 1 ,,p. That is “compression” refers to 
dimensionality reduction and not quantization which outputs bits for digital transmission. This 
is similar to the way considered, in particular, in lUTIl . 

B. Known techniques 

It is known that in the nondisrtibuted setting (in the other words, in the case of a single 
sensor only) the MSE optimal solution is provided by the Karhunen-Loeve transform (KLT) 
(131 , lfT4ll . (T8l , m. Nevertheless, the classical KLT cannot be applied to the above WSN since 
the entire data vector y = [yf,..., y^] T is not observed by each sensor. Therefore, several 
approaches to a determination of mathematical models for Q 1; ..., Q p and V have been pursued. 
In particular, in the information-theoretic context, distributed compression has been considered 
in the landmark works of Slepian and Wolf [fl9l , and Wyner and Ziv [[20tl . A transform-based 
approach to distributed compression and the subsequent signal recovery has been considered 
in urn, urn m, m, m, urn Intimate relations between these two perspectives have 
been shown in (251 . The methodology developed in |H71 . 11211 . (221 . (23l . fl24ll is based on 
the dimensionality reduction by linear projections. Such an approach has received considerable 
attention (see, for example, (5), 0, (H, ED, (26l, [[271). 

In this paper, we consider a further extension of the methodology studied in [1171 . (211, (221, 
(231 . (241 . In particular, in (171 . two approaches are considered. By the first approach, the fusion 
center model, V, is given in the form V = [P where V s is a ‘block’ of V, for 

j = 1 ,,p, and then the original MSE cost function is represented as a sum of p decoupled 
MSE cost functions. Then approximations to Qj and are found as solution to each of the 
p associated MSE minimization (MSEM) problems. We observe that the original MSE cost 
function and the sum of p decoupled MSE cost functions are not equivalent. This is because the 
covariance matrix cannot be represented in a block diagonal form in the way presented in IfTTl . 
Some more related explanations can be found in (28l . Therefore, the first approach in IfTTl leads 
to the corresponding increase in the associated error. The second approach in (171 generalizes 
the results in (231, (22l in the following way. The original MSE cost function is represented, by 
re-grouping its terms, in the form similar to that presented by a summand in the decoupled MSE 
cost function. Then the minimum is seeking for each VjQj , for j — 1,... ,p, while other terms 
V k Qk, for k = 1 ,..., j — 1 , j + 1 ,... ,p, are assumed to be fixed. The minimizer follows from 
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the known result given in If29l (Theorem 10.2.4). While the original MSEM problem is stated 
for a simultaneous determination of <2i,..., Q p and V, the approach in [H71 requires to solve p 
local MSEM problems which are not equivalent to the original problem. To combine solutions 
of those p local MSEM problems, approximations to each sensor model Qj are determined from 
an iterative procedure. Values of Qi,... ,Q P for the initial iteration are chosen randomly. 

The approach in |[24| is based on the ideas similar to those in the earlier references ifTTTl . [|2Tjl- 
GZL lf23lQ i.e. on the replacement of the original MSEM problem with the p+ 1 unconstrained 
MSEM problems for separate determination of approximations to Qj , for each j = 1..... jx and 
then an approximation to V. First, an approximation to each Qj , for j — 1,... ,p, is determined 
under assumption that other p—1 sensors are fixed. Then, on the basis of known approximations 
to Qi,, Q p , an approximation of V is determined as the optimal Wiener filter. Those p + 1 
problems considered in ll24l are not equivalent to the original problem. In [|24|. the involved 
signals are assumed to be zero-mean jointly Gaussian random vectors. Here, this restriction is 
not used. 

The work in [|2T|. |[22|. f[23ll can be considered as a particular case of |f24][. 

The method in m is applied to the problem which is an approximation of the original 
problem. It implies an increase in the associated error compared to the method applied to the 
original problem. Further, it is applicable under certain restrictions imposed on observations and 
associated covariance matrices. In particular, in lUTil . the observations should be presented in 
the special form y :i = 7/ ; x + v p for j = 1,... . p (where // y is a measurement matrix and v ; 
is noise), and the covariance matrix formed by the noise vector should be block-diagonal and 
invertible. It is not the case here. 


C. Differences from known methods. Novelty and Contribution 

The WSN models in lUTl . [H7I , ll2Tl . If22l . [[23il , [[24| are justified in terms of inverse matrices. 
It is known that in the cases when the matrices are close to singular this may lead to instability 
and significant increase in the associated error. Moreover, when the matrices are singular, the 
algorithms ifTTl . ll2T1 . Il22l . [[23H , ll24l may not be applicable. This observation is illustrated by 
Examples [2j [3] and [4] in Section VI below where the associated matrices are singular and the 


'in particular, it generalizes work HU to the case when the vectors of interest do not to be directly observed at the sensors. 
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method IfTTI is not applicable. Although in If24l . for the case when a matrix is singular, it is 
proposed to replace its inverse by the pseudo-inverse, such a simple replacement does not follow 
from the justification of the model provided in [f24j. As a result, a simple substitution of the 
pseudo-inverse matrices instead of the inverse matrices may lead to the numerical instability as 
it is shown, in particular, in Example [4] in Section VI In this regard, we also refer to references 
0, a, ea where the case of rank-constrained data compression in terms of the pseudo-inverse 
matrices is studied. 

Thus, methods in liTTfl . lETll . Il22l . |[23l . If24l . ifTO are justified only for full rank matrices used 
in the associated models. This is not the case here. On the basis of the methodology developed 
in 0, 0, [i30l , our method is rigorously justified for models with matrices of degenerate ranks 
(Sections |III-A| and |VIII[ ). 

Further, the second approach in ITF71 is, in fact, the block coordinate descent (BCD) method 
ED. The BCD method converges (to a local minimum) under the assumptions that the objective 
function and the space of optimization need to be convex or the minimum of the objective 
function is uniquely attended (see OTIl and lf32l . p. 267). Those conditions are not satisfied for 
our method (as for the method in lUTTl as well). Unlike the BCD method OTIl used in IfTTI . 
the maximum block improvement (MBI) method [HI, [O avoids the requirements of the BCD 
method. The MBI method guaranties convergence to a coordinate-wise minimum point which 
is a local minimum of the objective function. Therefore, the approach proposed in this paper is 
based, in particular, on the idea of the MBI method. 

As distinct from the technique in IfTTI our method is applied to the original minimization 
problem, not to an approximation of the original problem as in IfTTI . It allows us to avoid the 
increase in the associated error. We also do not impose any of the restrictions on our method as 
those in IfTTI mentioned in Section I-B above. In particular, we do not assume that the covariance 
matrix formed by the noise vector should be block-diagonal. 

The methods in IfTTI . lf2TI . If22l . If23l . Il24l have been developed under assumption that exact 
covariance matrices are known. A knowledge of exact covariance matrices might be a restrictive 
condition in some cases, in particular, when matrices are large. In this paper, these difficulties 
are mitigated to some extent. 

Key advantages of the proposed method are as follows. The method represents a combination 
of the generic Karhunen-Loeve transform (KLT) 0 (based on the work in 0) and the MBI 
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method [T], 0. Therefore, it unites the merits of both techniques 0, [[5! and 0]|, 0. As a result, 
our approach provides, in particular, the minimal MSE at each step of the version of the MBI 
method we use. The WSN model is represented in the form called the multi-compressor KLT- 
MBI transform and is based on the ideas which are different from those used in known methods 

[0, 0, [ill, m, 03, m, [123, 0D, [EH, [EH- The multi-compressor KLT-MBI is given 

in terms of pseudo-inverse matrices and, therefore, it is numerically stable and always exists. 
In other words, the proposed WSN model provides compression, de-noising and reconstruction 
of distributed signals for the cases when known methods either are not applicable (because 
of singularity of associated matrices) or produce larger associated errors. This observation is 


supported, in particular, by results of simulations in Section VI 


D. Notation 

Here, we provide some notation which is required to formalize the problem in the form 
presented in Section II-B below. Let us write (Q, £, //) for a probability spaceQ We denote by 
x G L 2 (fl, M r ") the signal of interest (a source signal to be estimated) represented as x = 
[x^,..., x( m )] T where x® G L 2 {f 2, M), for j = 1,..., m. Futher, yy G L 2 (Q, M ni ),..., y p G 
L 2 (Q, M” p ) are observations made by the sensors where ri\ + ... + n p = n. In this regard, we 
write 

y=[yf,....y„T and y= [y (1) ,....y'">] T (1) 

where y' k} G L 2 (Q, M), for k — 1,..., n. We would like to emphasize a difference between y ? 
and y (fc ); in ([!]), the observation y 3 , for j — 1,... ,p, is a ‘piece’ of random vector y (i.e. y 3 is a 
random vector itself), and y^ k \ for k — 1,..., n is an entry of y (i.e. y^ is a random variable)^] 
For j = 1,... ,p, let us define a sensor model Qj : —> L 2 (f2,M rj ) by the relation 


[Qj(yj))H = QjMu)] 


( 2 ) 


where Qj G 


r = ri + ... + r p , where r,- < n 


(3) 


= {a;} is the set of outcomes, E a < 7 —field of measurable subsets of and ji : E —> [0,1] an associated probability 
measure on E with /r(f2) = 1. 

’The space L 2 (Q,W m ) has to be used because of the norm introduced below in |5i. 

4 Therefore, y, = [y('*o+m+...+^_ 1 ) ) _ _ _ ^ y (m+...+nP]T where j = ^ p and nQ = L 
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Let us denote u j = Qfyf and u = [uf,..., uJ] T , where for j = vector u, e 

L 2 (Q,W j ) represents the compressed and filtered information vector transmitted by a jth sensor 
Qj to the fusion center V. A fusion center model is defined by V : L 2 (Q,W) —>■ L 2 (f2,M m ) so 
that 

[P(u)](o;) = P[u(w)], (4) 

where P e M mxr and u e L 2 (Q, M r ). To state the problem in the next section, we also denote 

E [ll x Hlli] : = ll x lln= [ ||x(w)|| 2 d/r(a;) < cx), (5) 

Jii 

where ||x(a-?) || 2 is the Euclidean norm of x(u) G M r ". For convenience, we will use notation 
||x||^, not E [||x(o;)|||], to denote the norm in (|5]). 


II. Formalization and Statements of the Problems 
A. Preliminaries and Formalization of the Problem 

For the WSN depicted in Fig. |T] the problem can be stated as follows: Find models of the 
sensors, Q 1} ..., Q p , and a model of the fusion center, V. that provide 


mm 

V,Qi,-.,Q p 


X-P 


Qi(yi 


Qp(yp 


( 6 ) 


under the assumption that Qi,... ,Q P and V are given by ([2]) and ([4]), respectively. The model 
of the fusion center, V, can be represented as V = [P i,...,P p ] where Vj : Ififl. M T ' J ) — y 
L 2 (Q,W n ), for j = 1,... ,p. Fet us write 

Qi(yi) 


mm 

Qi,...,Q P 


X- 


Q p (yp 


n 


= ||x - [PiQi(yi) + ... + V P Q P ( y p )]|| n 

Qi,—,Qp 


(7) 


where y = [yf,..., yW. 
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B. Statement of the Problem 

For j — 1,... ,p, let us denote F, = VjQj. We also write F = F\ ,..., F v ]. Then the WSN 
model can be represented as 

p 

F(y) = VQ{y) and F( y) = (8) 

3= 1 


Denote by 7 Z(m, n, k) the variety of all m x n linear operators of rank at most k. For the sake 
of simplicity we sometimes will also write 7 Zk instead of IZim. n. k). The problem in ([7]) can 


equivalently be reformulated as follows: Find Fi G Pirn, n,, rq), 
solve 


mm 

tF\ £'R-r 1 ,•••, tFp&Zr p 


3 =1 


F p G 7 Z(m, n p , r p ) that 

(9) 


Recall that r 3 < rij, for j — 1,... ,p. 

Further, to simplify the notation we will use the same symbol to denote an operator and the 
associated matrix. For example, we write Qj to denote both the operator Q ; and matrix Qj 
introduced in ([2]). Similarly, we write Pj to denote both operator Vj and matrix Pj G 
introduced above, etc. In particular, by this reason, in 0 we will write Fj , not F 3 . 


C. Assumptions 

For x represented by x = [x^\ ..., x^ m )] T where x^ g L 2 (f2, M) we write 

£[xy T ] = E xy = 6 K rox ”, (10) 

where (x (j) y (fc) ) = / ^\uj)y^ k \uj)dp{uj) and y = [y (1 \ ... ,y^] r . 

Jn 

The assumption used in the known methods []6]|, 0, flU, iTTOll . iflTI . lUTll . [I2T1 . Il22ll . Il23l . 
(241, ll26l - (271 is that covariance matrices E xy and E yy are known. At the same time, in many 
cases, it is difficult to know exact values of E xy and E yy . For instance, if it is assumed that 
signal y is a Gaussian random vector then the associated parameters p and a 2 for matrix E yy 
are still unknown (see (24l as an example). Therefore, we are mainly concerned with the case 
when estimates of E xy and E yy can be obtained. Methods of estimation of matrices E xy and 
E yy were studied in a number of papers (see, for example, (33l . (34l . (35l . (361 . (37l . (381 . 
(39l . (40l ) and it is not a subject of our work. In particular, samples of training signals taken 
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for some different random outcomes c o might be available. Some knowledge of the covariances 
can also come from specific data models. Then the aforementioned covariance matrices can be 
estimated. 


In the following Section III we provide solutions in terms of both matrices E xy , E yy and 


their estimates. The associated error analysis is given in Sections [V] and VIII 


D. Solution of Problem ([9]) for Single Sensor, i.e. for p = 1 

Here, we recall the known result for the case of a single sensor, i.e. for p = 1 in ([6]) and ([9]), 
and provide some related explanations which will be extended in the sections that follow. 

The Moore-Penrose generalized inverse for a matrix M is denoted by Mf We set M 1//2 t = 
(M 1 / 2 ) t , where M 1//2 is a square root of M, i.e. M = M 1//2 M 1//2 . Then the known solution (see, 
for example 01) of the particular case of problem (|9]), for p = 1, is given by 


F * = K,^.) 172 ],, ( B J*) 1/2 + Mi [/ - Eis,(Ei yi y/*] 


(ii) 


where symbol [-] ri denotes a truncated singular value decomposition (SVD) taken with ry first 


nonzero singular values and matrix M 1 is arbitrary. The expression (11) represents the Karhunen- 
Loeve transform as given, for example, in [0. 

In ( [IT] ), in particular, Mi = O where O is the zero matrix. Then ()\ and l\ that solve ([6]) 
follow from a decomposition of matrix [Ex^fE^J 1 ^ 2 ] r (E^ yiyi ) 1 / 2 , based on the truncated 
SVD, in a product P\Q\ of ry x n matrix Q i and m x r\ matrix Pi , respectively. 

We will extend this argument in Section HI below for a general case with more than one 
sensor. 


III. Main Results 

A. Greedy Approach to Solution of Problem ([9]), for p = 1,2,... 

We wish to find Fi,... ,F P that provide a solution to the problem ([9]) for an arbitrary finite 
number of sensors in the WSN model, i.e. for p — 1, 2,... in (|9]). To this end, we need some 


more preliminaries which are given in Sections III-A1 and III-A2 that follow. The method itself 


and an associated algorithm are then represented in Section III-A3 
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1) SVD, Orthogonal Projections and Matrix Approximation: Let the SVD of a matrix C 6 
I mxs be given by 


C = Uc^cVc, 


( 12 ) 


where Uc € K mxm and Vc G M sxs are unitary matrices, Ec = diag(cri(C),..., cr min ( mjS )(C')) e 
K raxs is a generalized diagonal matrix, with the singular values a\ (C) > a 2 (C) > ... 0 on the 
main diagonal. Let U c — [ui u 2 ■ ■ ■ u m \ and V c = [r : i v 2 ... v s ] be the representations of U 
and V in terms of their m and s columns, respectively. Let 

rank C rank C 

L c = UiUG RmXm and R c= J2 ViV i G RSXS (13) 

2—1 2=1 

be the orthogonal projections on the range of C and C T , correspondingly. Define 

r 

Cr = [C] r = °i(C)u iV J = Uc r ^c r Vl e R mxs (14) 

2=1 

for r = 1,..., rank C, where 

U Cr = {uiu 2 ... Ur], Ecv = diag(ai(C),..., a r (C )) and V Cr = [nr v 2 ... v r ], (15) 


For r > rank C, we write C (r) — C(— C' ran k c)- For 1 < r < rank C, the matrix C (r) is 
uniquely defined if and only if a r (C) > cr r+ \ (C). 

Consider the problem: Given matrices S :j and G t , for j — 1,... ,p, find matrix Fj that solves 


min \\Sj — FjGj \\ 2 , for j 

r j tMrj 




(16) 


The solution is given by Theorem |T] (which is a particular case of the results in 0, fiSl) as 
follows. 

Theorem 1: Let Kj = Mj (/ — Lq :i ) where Mj is an arbitrary matrix. The matrix Fj given 


by 


F) = [S,Rc,] ri G](I + Kj), for j 




(17) 


is a minimizing matrix for the minimal problem (16). Here, Rg : and [-] rj are defined similarly 
to ([13]) and (JT4]), respectively. Any minimizing Fj has the above form if and only if either 


Tj > rank (. SjR Gj ) 


(18) 


August 20, 2015 


DRAFT 



11 


or 


1 < rj < rank ( SjR Gj ) and a rj {SjR Gj ) > a r+1 (SjR Gj ), 


(19) 


where a rj (S ) R ( ; i ) is a singular value in the SVD for matrix SjR Gy 

Proof: The proof follows from 01, [j5]. ■ ■ 

We note that the arbitrary matrix Mj implies non-uniqueness of solution ( [T7| ) of problem ( [T6] ). 
2) Reduction of Problem ([9]) to Equivalent Form: We denote F — [F 1; ..., F p ], where F E 
R mxn and Fj G W nXnj , for all j = 1.... .p, and write || ■ || for the Frobenius norm. Then 

2 

= ||x - F(y)||n = tr {E xx - E xy F T - FE yx + FE yy F T } 


x -5^( yj 

3 = 1 


- WE^E^n 2 + ||- FE'*\\ 


( 20 ) 


Denote by M(m, n, k ) the variety of all m x n matrices of rank at most k. For the sake of 
simplicity we also write k). 


In (20), only the last term depend on F u ..., F p . Therefore, (|9j) and (20) imply 


mm 

J r i£'R-r 1 r --,J r p£'R-r p 


3 =1 
-, 1/2 


mm 


-FlGlRri ,...,-FpEM rp 


Let us now denote H — E, v ( El -/)' and represent matrix E^ 2 in blocks, Ef/ — [Of , ... G 


(2D 


yy 


nl/2 _ r^ T 


where G :) G M n ' xn , for j = 1,... ,p. Then in (21), 


(p] T 


- FEl ' 2 = H 

3 = 1 


Therefore, (20) and (|22|) imply 


mm 

J r i£'IZr 1 ,...,J r p£'R. rp 


p 

2 

p 

x -i>w 

i=i 

= min 

ft 

H~Y _; TjGi 

i=i 


On the basis of (23), problem ([9]) and the problem 

h - £ F = G > 


mm 

FiGM ri 


3 = 1 


( 22 ) 


(23) 


(24) 


are equivalent. Therefore, below we consider problem (24). In (24), for j — 1,... ,p, we use the 
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representation 


H-J2 FjG, 

3 = 1 


— \\Sj — FjGj \\ 2 , 


(25) 


where Sj = H — ^ F ?; Gj. 


1=1 
*7 G 


3) Greedy Method for Solution of Problem (24): A solution of problem (24) is based on the 
idea of the MBI method [[0, fill which is a greedy approach to solving optimization problems. 


The advantages of the MBI method have been mentioned in Section I-C To begin with, let us 
denote 

2 


Vl,...,r p — 


x ... x M rp , F — (Fi ,..., F p ) e 




and {(F) = 


H-Y. F i G > 

3 = 1 


The method we consider consists of the following steps. 

1st step. Given F (0) = (F^,..., F^) e compute, for j = 1,... ,p, 

p 

. L) = G\ (T 4- K A wherp — H — ' 


G](I + K, ), where sf = H - ^ F t (0) G', : . 


Z=1 


Note that F (1) is the solution of problem (16) represented by (17). A choice of i^ 0) 


is considered in Section [TV] below. 
2nd step. Denote 


p(i)_ ( p(0) p(0) p(l) p(0) 


f(0) 

I - 1 p 


select F, 1 such that 


/(A 1 ’) = _ min {/(ff),..../(A 1 ’). ■ ■ ■, /(A 11 : 

F i •■■■Tp j 


and write 


(26) 


^(°) 

i ■ ■ ■ i r p ) 


(27) 


(28) 


^ (1) = F ( k\ 


where we denote F (1) = (F^\ FjP) G 




(29) 


Then we repeat procedure (26)-(29) with the replacement of F <(,) by F (1) as follows: Given 
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F (1) = (F^\ Fp 1 ^), compute 


p j 2 ) — 


S^’Rc, 


G](I + Kj), where sf } = H - F^Gi, 


1=1 

*7 


( 2 ) 

select F k that satisfies (28) where superscript (1) is replaced with superscript (2), and set 

F { 2) = F ( k \ 

This process is continued up to the gth step when a given tolerance e > 0 is achieved in the 
sense 


| /(F ( «+ 1) )- /(F ( »))| < e> for, = 1,2,.... 


(30) 


It is 

summarized as 

follows. 






Algorithm 1: Greedy solution of problem (24 

) 

Initialization: F {() \ 

H, S u .. 

S'p and e 

> 0. 




1 . 

for q = 0,1, 2, 







2. 

for j = 1 , 2, 

-,P 






3. 

S ( f ] = H 

- E F i’ >G ‘ 







1=1 
*7 H 






4. 

p(<z+i) _ 

3 

q (?) d 

3 

Vi 





5. 

p{q+l) _ 

(ff”,.... 

jq(l) £3(9+1) rril) p(9) 



6. 

end 







7. 

Choose F ( ' g+1 ' > 

= F^ +1) where F 

(9+1) = 


..,F P (,?+1) ) and 

F^ +1) is such that 


f(F* +1) ) = min {/(ff +I >),... 

r 1 

,/(A ,+1> )... 


7. 

If |/(F (9+1) ) - 

- f(F^) 

1 <e 





8. 

Stop 







9. 

end 







10 

end 








Algorithm 1 converges to a coordinate-wise minimum point of objective function f(F) = 
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p 

2 


- £ Ffi, 

3=1 

which is its local minimum. Section 

VIII-A 


provides more associated details. 


Remark 1: For p — 1, the formulas for Fj in ([17]) and (26) coincide, and take the form 


Fi = [HR Gl 1 G\(I + K x 


(31) 


where G 1 = E^, K 1 = M 1 (/ — L Gl ), L Gl = an d HR Gl = E xyi {E\^fl\ 

Thus, Fi in ( |3T[ ) coincides with Fi in ( [IT] ) (since Mi is arbitrary). In other words, the KLT 
represented by ( |TT| ) is a particular case of the expressions in (17) and (26). For this reason, ( [17] ) 
and ([26]) can be regarded as extensions of the KLT to the case under consideration. Therefore, 


the transform represented by F ; (y ; ) where F\..... F p solve (16) can be regarded the multi- 


3 = 1 


compressor Karhunen-Loeve-like transform (or the multi-compressor KLT). Further, Algorithm 
1 represents the version of the MBI method that uses the multi-compressor KLT. Therefore, the 

p 

WSN model in the form F • 7 ' 11 


(y ) where F 1 ( ' 9+1 \ ..., Fp 111 j are determined by Algorithm 1 


3 = 1 


can be interpreted as a transform as well. We call this transform the multi-compressor KLT-MBI. 

Remark 2: In practice, an exact representation of matrices E xy and E yy might be unknown. 
Their estimates, E xy and E yy , can be obtained by known methods Il33l . [[34ft . Q5I . 0511 . 071 . 
], m. In this regard, we denote H = E xy (El y 2 y and F,^ 2 = [Gj\... ,Gj] T where 

j2 

, for j — 1,... ,p, is a block of E yy . Then in (26), Gj, H and Gi should be replaced 
with G y , II and G t , respectively. In this case, steps (26)-(30) of the method and Algorithm 
1 are the same as before but F < ' 9+1 ^ = (F^ 9+1 \..., Fp Q+1 ^) should be denoted by F 


Gj G R niXn 


(9+1) 


^^7(9+1) ^( 9 +lb 


B. Models of Sensors and Fusion Center 

The models of sensors Qi,...,Q p and the fusion center P = [P 1 ,...,P p ] of the WSN in 
Fig. jT| follow from the formula Fj = PjQj introduced in ([8]). By the proposed method, Fj is 
represented by Fj 1 ' given by Algorithm 1 above. Therefore, the mathematical model of the 
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WSN is given by 


X>j ,+1) (y/) = E^ ,,+1| Q5 ,+1| (yj) = p 1s+1) 


3 = 1 


3 = 1 


' 0S 8+1) (yi) N 
v Q? + 1 , (y P ) , 


(32) 


Here, Fj 9+1) = Pj 9+1) Q [q+1) and P^ +1 ) = [P 1 (9+1) ,..., P p (9+1) ], Then qJ 9+1) £ M rj ' Xnj ‘ and 
pjy+!) e f™xr i f 0 u 0W from the SVD decomposition of matrix F^ l+V> taken with r 3 first nonzero 
singular values. This procedure has been described in Section [n-D| for the case of only one sensor, 
i.e. for p = 1. 

For the case when instead of matrices E xy and E yy their estimates E xy and E yy are used, the 
models are constructed by the similar procedure with the replacement of Pj 9+1 ^ by F J - 9+1 (see 
Remark [2] above). 

Note that the proposed WSN model also provides de-noising of observations yi,... ,y p . 


IV. Determination of Initial Iterations 

To start Algorithm 1, the values of initial iterations pj 0) and Fj should be defined. It is 
done as follows. Let us denote x = [xf,... , xJ] T where x* £ L 2 (f2,M mi ), i = 1 ,p and 
mi + ... + Trip = m. Suppose that matrix P £ M mxr is given by P = diag(Pn,..., P pp ) where 
Pjj £ M. mjXrj , for j = 1,... ,p. Then 


Xl 


diag(Pn,.. .,V PP ) 


Qi(yi) 


Qp( y F 


xi -Pi(yi) 


x P - P P (y P ) 


E 

3 = 1 


l x i -^-(yjOlln 


(33) 


and 


v 


mm 


l x i - ^i(y,-)||n = min ||xj - T y 


T\&n r ~ 1 ~..~T p &n rv ^ "~~ J JX 9 ^ F j eTi~ rj 


jj wa¬ 


rn 


3 = 1 


3 = 1 


August 20, 2015 


DRAFT 




















16 


As a result, in this case, problem ([9]) is reduced to the problem of finding T t that solves 


A0 r IN-^WIIn 

J 3 


(35) 


for j — 1,... ,p. Its solution is given by ( [IT] ) in Section |II-D| above, i.e. by 


Fj = 


(K in y /2 + W - 


(36) 


where is an arbitrary matrix. Then the initial iterations for Algorithm 1 are defined by 

F- 0) = Fj, for j = 1,... ,p. 


(37) 


Similarly, when instead of matrices E xy and E yy their estimates E xy = {E x . y . }f J=1 and E yy = 
{E yiyj } p ij=l are used, the initial iterations are defined by 


p )°) — 


^WAV 172 ! (4„) i/2 +m - ei'keip'% 


VjVi 


y.:>h ' VjVj ■ 


(38) 


where E x . y and E y . y . are blocks of E xy and E yy , respectively. 


"ViVi 


V. Error Analysis: A Posteriori Associated Errors 

For F (9+1) = [F-f 

proposed WSN model is represented as 


..., Fp 1 n determined by Algorithm 1, the error associated with the 


I|X - [F, (,+1) ,.... ^ + 1 ) ](y)lln = ll^ll 2 - l|B,(^ 2 ) t l | 2 + I \E v (Elf)l - 

where F^ q+l) = [F 1 (,+1) ,..., F p (9+1) ], 

For F^ l+V) described in Remark 2 the associated error is given by the similar expression: 

||x-[F 1 ( ’ +1) ,...,F'»+ 1 '](y)||? ! = ||E^ 2 || 2 -||F !TO (^ 2 ) t || 2 + ||F I „(F ; ;' 2 ) , -F ( «+ 1 )E;/ 2 || 2 , (39) 


where = [i^ 9+1 \..., Fp 1+l> ], The above formula (39) is used in our simulations repre¬ 

sented in the following section. 


VI. Simulations 

Here, we wish to illustrate the advantages of the proposed methodology with numerical 
examples carried out under the assumption that either covariance matrices E xy , E yy or their 
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estimates are known. The assumption that only the covariance matrices are known is similar 
to that used in (6), 0, [®, OH, CTD, El, GUI, [E3, fl23j, GiU, (261, El- In particular, the 
estimates can be obtained from samples of training signals. In many situations, the number of 
samples, s, is often smaller than the dimensions of the signals x and y, which are m and n, 
respectively l!38l . At the same time, it is known that as s —> oo, the ergodic theorem asserts 
that the estimates converge to the true matrix values ll39l . PTO1 . In particular, for large s, the 
estimates of the covariance matrix have been considered in [134-1 . [f35l . It is interesting to compare 
our simulation results for the cases when s is ‘relatively’ small and ‘relatively’ large. In the 
examples that follow, both case are considered. 

A comparison with known methods [fill . ETl . El . El . El . El is a s follows. The method 
El represents a generalization of methods El . El . El and therefore, we provide a numerical 
comparison with method El which includes, in fact, a comparison with methods El . El- 
El as well. Further, covariance matrices used in the simulations associated with Figs. [3] (a), (c) 
and Figs. [4] (b), (c), (d) are singular, and therefore, method ifTTl is not applicable (in this regard, 
see also Section |II-C| ). Therefore, in Figs. [3] (a), (c) and Figs. [4] (b), (c), (d), results related 
to algorithm in IfTTl are not given. By the same reason, the method presented in IfTTl is not 
applicable as well. Moreover, the method in IfTTl is restricted to the case when the covariance 
matrix formed by the noise vector is block diagonal which is not the case here. 

In the examples below, different types of noisy observed signals and different compression 
ratios are considered. In all examples, our method provides the better associated accuracy than 
that for the methods in IfTTl . El (and methods in El . El . El as well, because they follow 
from El ). 

Example 1: We start with an example similar to that considered in El assuming that a WSN 
has two sensors and the observations y 1 and y 2 are represented by 


yi = x + fi, and y 2 = x + £ 2 , 


(40) 


where x e L 2 (fl,M 3 ), £] e L 2 (fl,M 3 ) and £ 2 6 L 2 (fl,M 3 ) are Gaussian independent random 

0.585 0.270 0.390 
0.270 0.405 0.180 


vectors with the zero mean. Let E TT = 


0.390 0.180 0.260 

where ay = 0.2 and a 2 = 0.4, and J 3 is the 3 x 3 identity matrix. Then E xy = [E xx E xx ] and 


and E{^ = cr 2 / 3 , for j = 1, 2, 
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Eyy ~ 


7j xx + I 3 


E x 


E X x E xx + 0"2^3 
achieve tolerance e = 0.146. The ac 


. For ri = r 2 = 1, Algorithm 1 requires three iterations to 

lievable tolerance of methods ifTTI . Il24l . for r 1 = r 2 = 1, is 
18% worse, e = 0.173, and it is not improved after the initial iteration proposed in 1(241 . 
Example 2: Let us consider the case of a WSN with two sensors again where, as before, 


yi = x + £i, and y 2 = x + £ 2 , (41) 

where x G L 2 (f2,M m ), G L 2 (f2,R m ) and £ 2 G L 2 (f2,R m ), and y : and y 2 are noisy versions 
of the source. Unlike Example jT] we now assume that covariance matrices E xy and E yy are 
unknown. Therefore, their estimates E xy and E yy = j E yiVj |, i,j = 1,2, should be used. To 
this end, estimates E xy and E yy have been determined from the samples of training signals as 
follows: 

E xy = l -[XY^XYj] and E wy . = -Y^ lortj = 1,2. (42) 

Here, X G M mxs has uniformly distributed random entries and, for i = 1, 2, 

Y^X + nTi, (43) 

where cq G M and T, G WL'"' XS has random entries, chosen from a normal distribution with mean 
zero and variance one. Diagrams of typical errors associated with the proposed Algorithm 1 and 
known methods IfTTI . Il24ll are given in Figs. [3] (a), (b). 

Note that in Figs. [3] (a), (b), the obtained results are illustrated for different compression ratios 
Cj = r j/rij where j — 1,.. . ,p. The compression ratios Cj = 1/5, for j = 1,2, 3, used to obtain 
the results represented in Fig. [3] (b) are smaller than those in Fig. [3] (a), c\ =3/5 and c 2 = 7/10. 
This is a reason for the error magnitudes represented in Fig. [3] (a) being smaller than those in 
Fig. [3] (b). This observation also holds for other examples that follow. Further, in Fig. [3] (b), 
due to large sample size, s = 10000, the estimate of matrix E yiVj is very close to its true value 
which is the identity. By this reason, iterations of our method are similar to each other and the 
associated errors are similar for almost all iterations. The same effect holds for methods IfTTI . 


Example 3: Here, we consider the case when observations y, and y 2 are very noisy, i.e. 
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reference signal x is significantly suppressed. To this end, we do not assume that Y t e M mxs 
is represented in the form ([43]) but it has random entries, chosen from a normal distribution 


with mean zero and variance one. We also use E xy and E yy = | E yiyj j, i,j = 1,2 in the 

is as in the above Example [2j For m = n\ = N 2 = 
0.086 0.439 0.857 0.904 
0.074 0.574 0.386 0.429 
0.4660 -0.1260 0.3870 0.3290 


form (42) as before where X e 


2 and s = 4, examples of those matrices are X = 


, = 


0.284 -0.942 0.067 0.222 

-2.206 0.514 -1.293 -0.686 


and Yo = 


0.6880 -0.4690 -0.9420 -0.5630 


For the case of two sensors (i.e. for p — 2) and for the above samples X, Y\ and Y 2 , the 
errors associated with the proposed method and known method |j[24l are given in Fig. [3] (c). 

For larger magnitudes of m,ni, n 2 , n 3 , s and ry, for j = 1, 2, 3, the errors associated with the 
proposed method and known methods are represented, for the case of two and three sensors (i.e. 
for p = 2 and p — 3, respectively), in Fig. [3] (d) and Fig. [4] (a). 

Example 4: In this example, we consider the case when observations are corrupted by noise 
in the way which is different from those in Examples [2] and [3] Namely, we assume that, for 
3 = 

Yj = A'X + si (44) 


where Aj : L 2 (Q. R r ") ->■ L 2 (Q. R m ) is a linear operator defined by matrix Aj e R mxm with 
uniformly distributed random entries, and ^ is a random noise. Samples of x and £,■ are simulated 
as matrices X e M" IXS and Tj e M mx ' s , respectively, where crj e M, such that X has uniformly 
distributed random entries and Tj has random entries, chosen from a normal distribution with 
mean zero and variance one. 

The errors associated with the proposed method and the known method, for the case of two 
and three sensors (i.e. for p = 2 and p — 3, respectively), and different choices of m, n 3 . s and 
Tj, for j — 1,2 and j = 1, 2, 3, are represented in Figs. [4] (b), (c) and (d). 

Note that method [|24| is not numerically stable in these simulations. We believe this is because 
of the reason mentioned in Section II-C1 

Example 5: In the above examples, we used estimates of training signals, not training signals 
themselves. Here, we wish to illustrate the obtained theoretical results in a different way, by a 
comparison of a training reference signal with its estimates obtained by our method and known 
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methods. To this end, we simulate the training reference signal x by its realizations, i.e. by a 
matrix X £ M mxk where each column represents a realization of the signal. A sample A" £ R mx ’ s 
with s < k is formed from X by choosing the even columns. To represent the obtained results 
in a visible way, signal X £ M' mxfc i s chosen as the known image Lena given by the 128 x 128 
matrix - see Fig. [ 5 ] (a), i.e. with m, k = 128. Then X £ M 128x64 . 

Further, we consider the WSN with two sensors, i.e. with p — 2, where the observed signal 
Yj, for j = 1,2, is simulated as follows: 


Yj = Aj*X + d T 


J d J 


where A~ £ 


3128x128 


has uniformly distributed random entries, T 3 £ 


B 128 xl 28 


has random 


entries, chosen from a normal distribution with mean zero and variance one, Aj*X represents 
the Hadamard matrix product, and a\ = 0.2 and cr 2 = 0.1. Estimates E xy and E y . y . are used in 
the form (42) where sample Yj £ M 128x64 is formed from Y l by choosing the even columns. 


For r 1 = r 2 = 64, the simulation results are represented in Figs. [5] and [6j Our method and 
known methods in ll24ll and iflTl ) have been applied to the above signals with 50 iterations each. 
The associated errors are evaluated in the form \X — X\ where X is the reconstruction of X 
by the method we use (i.e. by our method or methods in IfTTII and ||24|0 . 

Similar to the other examples, Figs. [5] and [6] demonstrate a more accurate signal reconstruction 
associated with the proposed method than that associated with known methods. 


VII. Conclusion 

We have addressed the problem of estimating an unknown random vector source when the 
vector cannot be observed centrally. In this scenario, typical of wireless sensor networks (WSNs), 
distributed sensors are aimed to filter and compress noisy observed vector, and then the com¬ 
pressed signals are transmitted to the fusion center that decompress the signals in such a way 
that the original vector is estimated within a prescribed accuracy. The key problem is to find 
models of the sensors and the fusion center in the best possible way. 

We proposed and justified the method for the determination of the models based on a combina¬ 
tion of the solution [O of the rank constrained least squares minimizing problem (represented by 
([9]) in Section II-B) and the maximum block improvement (MBI) method [[0, [12 ■ The proposed 
method is based on the following steps. First, we have shown how the original problem can be 


August 20, 2015 


DRAFT 




21 


reduced to the form ( f2T[ ) (Section |III-A2[ ) that allowed us to use the approaches developed in 
id, a. As a result, under the assumption that the associated covariance matrices or their 
estimates are known (from testing experiments, for example), the procedure for determining 
models of the sensors and the fusion center is given by Algorithm 1 (Section III-A3[ ). 

The obtained optimal WSN model represents an extension of the Karhunen-Loeve transform 
(KLT) and has been called the multi-compressor KLT-MBI. The known KLT follows from the 
multi-compressor KLT-MBI as a particular case. The models of the sensors and the fusion center 
have been determined in terms of the pseudo-inverse matrices. Therefore, the proposed models 
are always well determined and numerically stable. In other words, the proposed WSN models 
provide compression, de-noising and reconstruction of distributed signals for the cases when 
known methods either are not applicable or produce larger associated errors. As a result, this 
approach mitigates to some extent the difficulties associated with the existing techniques. Since a 
‘good’ choice of the initial iteration gives reduced errors, the special method for the determination 
of the initial iterations has been considered. 

The error analysis of the proposed method has been provided. 

Finally, the advantages of the proposed method have been illustrated with numerical experi¬ 
ments carried out on the basis of simulations with estimates of the covariance matrices. It has 
been shown, in particular, that the errors associated with the proposed technique are smaller than 
those associated with the existing methods. This is because of the special features of our method 
described above. 


VIII. Appendix 


A. Convergence 


Convergence of the method presented in Section III-A3 can be shown on the basis of the 
results presented in 0, [[2l as follows. 

We call F = (F 1 ,..., F p ) e a point in the space M ri) _ >rp . For every point F 6 M ri . 

define a set 


= {(F h Fj_ i)} x K r x {(F j+ 1 ,..., F p )} , for j — 1,... j p. 


A coordinate-wise minimum point of the procedure represented by Algorithm 1 is denoted by 
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F'=(F;,...,F;) wherfl 


F € («rg min f(F;,...,F;_ 1 ,F i ,F; +1 ,...,F;) \ . 


This point is a local minimum of objective function in (24), f(F ) = 


that F l ' g+1 ' ) in Algorithm 1 and F* defined by (45) are, of course, different. 
For F q> defined by Algorithm 1, denote 


j= 1 


(45) 


Note 


F = lim 

q—>oo 


(46) 


Note that because of ([H]), the sequence {F^} is bounded. 


Theorem 2: Point F defined by (46) is the coordinate-wise minimum of Algorithm 1. 

Proof: For each fixed F = (F\,.... F p ), a so-called best response matrix to matrix Fj is 
denoted by where 


e ) ar S f(Fi,-Fj-i,Fj,F j+ i,...,F p ) 

{ Pj£Rrj 

Let {F^} be a sequence generated by Algorithm 1, where F (q> = {F[ q \ ..., Fp l] ). Since each R r . 
is closed [41, p. 304], there is a subsequence {.F- 9< 4} such that (Ff 9s \ ..., Fp 9s - > ) —y (F *,..., F*) = 
F* as s —* oo. Then, for any j = 1, ...,p, we have 

/(FA •••, FA F". FA FA a /(FA •••, FA *f"FA FA 

\ + p(9s+l) p(<Js + l) p(?3 + l) zp(g s +l)\ 

— j Fl , iJFj I r j +1 1"-I r p ) 

f ( Zphla+l) eF9s+i) ri(9a+l) 7p(<Zs+l) zp(g s -|-l) N 

— J \ F 1 1 ) "j 1 r j+ 1 / 

By continuity, when s —» oo, 


/(J 7 T.....r*-i. F .F+i.....F) > /(F*. •••■ F-i. AF+i> -.F). 

which implies that above should hold as an equality, since the inequality is true by the definition 


The RHS in 1 45 i is a set since the solution of problem min f{F*,...,Fj-uFj,Fj +1 ,...,Fp) is not unique. 


6 There could be other local minimums defined differently from that in 1 45 .. 
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of the best response matrix <P ■ . Thus, F* is such as in (45), i.e. F* is a solution of the problem 


mm f(F;,...,F;_ 1 ,F j ,F' +1 ,...,F;), Vj = 1 


Fj&K rj 


Remark 3: Theorem 


still holds if the objective function is defined as f(F) = 


where F = ( F\ ,..., F p ), and H and G t are defined by Remark 2 In this case, the coordinate- 


H-Y.FG, 

3 =1 


wise minimum point is defined similar to that in (45) where symbols F*, F* and Fj should be 


~ ( q ) 

replaced with F , F* and Fj, respectively. More precisely, if F denotes the gth iteration of 
Algorithm 1 as described in Remark [2j then the following is true. 

Corollary 1: Point F defined by 


F = lim F 

q—toc 


(9) 


is the coordinate-wise minimum of Algorithm 1 in the case when covariance matrices E xy and 
E yy are replaced with their estimates E xy and E yy , respectively. 
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Iterations 

(a) Example [2] p = 2, m = rij = 10, s = 20, rj = 5 + j, 
<jj = 0.2 — 0.1 j, for j = 1, 2,. 



Iterations 

(c) Example[3] p = 2, m = rij = 2, s = 4, rj = 1, cn, = 1, 
fori = 1,2. 



Iterations 

(b) Example [2] p = 3, to = rij = 100, s = 10000, rj = 
20, (jj = 1, for j = 1, 2, 3. 



Iterations 

(d) Example [3] p = 2, m = rij = 100, s = 250, rj = 20, 
<jj = 1 , for i = 1 , 2. 


Fig. 3. Diagrams of MSE’s associated with the proposed method and the known methods versus number of iterations, for 
p = 2 (i.e. for two sensors) and p = 3 (i.e. for three sensors), and different choices of signal dimensions, m, rij, rj, sample 
size s and noise ‘level’ <jj. 
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Iterations 

(c) Example [4] p = 3, m = rij = 100, s = 400, rj = 25, 
o’.) = 0.1 j, for j = 1,2,3. 



Fig. 4. Diagrams of MSE’s associated with the proposed method and the known methods versus number of iterations, for 
p = 2 (i.e. for two sensors) and p = 3 (i.e. for three sensors), and different choices of signal dimensions, m, n,, rj, sample 
size s and noise ‘level’ <jj. 
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(c) Observed signal > 2 . 


(d) Estimate of A' by our method. 



20 40 60 80 100 120 20 40 60 80 100 120 

(e) Estimate of X by method CD (f) Estimate of X by method m 


Fig. 5. Illustration to Example [5] 
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(a) Error associated with our method. 



(b) Error associated with method ED- 



Fig. 6. Illustration to Example [5] 


(c) Error associated with method Hi- 
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